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1. INTRODUCTION 

. The behavior of fault-tolerant computing systems can be modeled 
as continuous- time Markov processes on large state spaces. Calculation 
of reliability is equivalent to computing transient probabilities of 
states of the Markov process corresponding to system failure. The state 
spaces are often very large, and thus efficient computational methods 
are required in order to calculate state probabilities. The CARE III 
approach has been developed to solve this problem; it is presented by 
Stiffler, Bryant, and Guccione [8] and further discussed by Trivedi and 
Geist [9]. The "randomization" technique is an alternate approach which 
is of considr 'ble interest in its own right and which will be useful in 
va3.idating the CARE III approach for systems with moderate state spaces. 

The randomization modeling and computational technique will be 
illustrated on a simplified model of a fault-tolerant system consisting 
of three components similar to one presented by Trivodi and Geist [9]. 
Figure 1 shows the behavior of a single component. Initially the component 
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The randomization technique for computing transient probabilities 
of Markov processes is presented. The technique is applied to a Markov 
process model of a simplified fault-tolerant computer system for illus- 
trative purposes. It is applicable to much larger and more complex 
models. Transient state probabilities are computed, from which reli- 
abilities are derived. A new accelerated version of the randomization 
algorithm is developed which exploits "stiffness" of the models to gain 
increased efficiency. A great advantage of the randomizaton approach 
is that it easily allows probabilities and reliabilities to be computed 
to any predetermined accuracy . 
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0 ; component OK 
A : active fault 

B : benign fault 

D : detected (and reconfigured) 

E : error 

F ; failure (error propagated) 

Figure 1. — Single component reliability 
model: state space, 

transitions, and rates. 


is fault-free, but after an exponential holding time with rate X an 
"active" fault occurs. Prom the active state the fault may become "be- 
nign" and later become active and continue alternating between active 
and benign. From the active state the fault may be "detected" (by diag- 
nostics) or generate an "error." This error may lead to detection of 
the fault and system reconfiguration or to "failure" of the system. 
Figure 1 shows the six states of a component, the possible transitions, 
and the rates at which they occur. We shall apply the randomization 
procedure to a system consisting of three independent components. The 
state space of the three-component system is shown in Figure 2; also 
shown are the possible transitions of the Markov process and their rates 
The model has 18 states and 31 transitions. The set {F, AA, AE, DF, FD 
XXA} is defined as "system failure," and the goal is to compute the 
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Figure 2. — A three-eoinponent system reliability model: 
state space, transitions, and rates. 

probability of absorption of the process in this set at time t^^j , the 
mission completion time. Note that this model is similar to one pre- 
sented by Trivedi and Geist f9, p. <^(6]. It has a small state space, but 
it serves well as an illustrative example. The randomiaation technique 
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can be efficiently applied to models with much larger state spaces. It 
should be useful in analyses of two fault-tolerant computing systems 
developed under NASA sponsorship: SIFT [10] and FXMP [5], 

The paper is structured as follows. Section 2 contains the stan- 
dard randomization algorithm for computing transient probabilities of 
Markov processes. In Section 3 a new accelerated version of the random- 
ization procedure is developed; it exploits the fact that models of many 
fault-tolerant computing systems are "stiff," i.e., the model has very 
fast and very slow transition rates. Section 4 gives computational re- 
sults for the standard and accelerated algorithms applied to ten differ- 
ent versions of the three-component model of Figure 2. Section 5 con- 
tains summary comments and a brief discussion of other approaches. Two 
appendices contain listings of FORTRAN programs of the two algorithms. 
Additional information on the randomization technique may be found in 
Gross and Miller [3]. 

2. THE STANDARD RANDOMIZATION ALGORITHM 

Let {x(t) , t ^ 0} be a continuous-time Markov process on a 
finite state space S = {l,2,..,,m} . The state probability vector at 
time t is denoted TT(t) = •••« ‘H’ (t)} , where Tr (t) = 

p([x(t) = s3 , s e S , Two different characterizations of the stochastic 
nature of (x(t) , t > 0} are useful: (i) the infinitesimal generator 

and (ii) a randomized Markov chain . 

All Markov processes can be characterized by an initial distribu- 
tion 2[(0) and an infinitesimal generator 
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and 




The transition rates which are depicted in Figures 1 and 

2. The infinitesimal generator Q seems to be the most natural way to 
describe the stochastic nature of the Markov models of fault-tolerant 
computing systems. 

Any Markov process on a finite state space can be represented as 
a discrete time Markov chain "randomized" by a Poisson process. Define 

A = max q.. (2.1) 

ieS ^ 


and 


P = Q/A + I , (2.2) 

where I is the identity matrix; P is a stochastic matrix. Let 
{Y^, n = 0,1,2,...} be a Markov chain on 5 with transition matrix P 
and initial distribution jr(0) . Let (N(t), t > 0} be a Poisson pro- 
cess witii rate A which is independent of (Y^^, n = 0,1,2,...; . Then 
{Y^(tj, t ^ O} is a Markov process witli generator Q and initial dis- 
tribution Tr(0) and hence is probabilistically identical to {x(t), t 5= 0 
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(The relationship between sample paths of {Y^, n = 0,1,2,..,} , 

{N(t), t > 0} , and shown in Figure 3.) This con- 

struction makes it possible to compute transient probabilities of a 
Markov process with generator Q from transient probabilities of a 
Markov chain Y with transition matrix P and a Poisson process N 
with rate A , The transient probabill •*ics of Y are denoted ^(n) = 
C<{)^(n), ’ where <l>gvn) = P(Y^ = s) , s e S . The 

randomization formula is 


pCx(t) = s} = I pCx(t) = s I N(t) = n} pCN(t) = n} 
n=0 


I P(Y^ = s) pCN<t) = n] 

n=0 " 


or equivalently, 


“ “At^.^.n 

aC:) - I w 

n=0 


(2.3) 


See Gross and Miller [3] for additional discussion and details. (Equa- 
tion (2.3) can also be found in finlar [1, p. 259].) 

The infinite series in Equation (2.3) must be truncated for com- 
putational purposes. Let 


T(e,t) = minfk: P([N(t) > k^ < e] = min k: ^ > 1 - 

■' n=0 "• 


(2.4) 


whf*re e equals the acceptable error (specified by the user). The com- 
put.i>:'ional version of Equation (2.3) is 


T(c,t) -At..,.n 

- I ♦(.-) ° ’ 

n=0 


(2.5) 


Truncation of the infinite series involves a probability loss of at most 
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n 


N(t) 


8 

7 

6 

5 

4 

3 

2 

1 

0 


X(t) = 
^N(t) 


Figure 3. — Example oi; randomization; Markov chain Y^, Poisson process 
N(t), and randomized chain X(t). 


e , thus all probabilities (oJi states or subsets o£ states) will have an 
error between -e and 0. Note that the randomization formula (2.5) 
reduces the calculation of transient probabilities of a Markov process 


_ 7 _ 


T-469 


to thoso of a Markov chain and underlying Polsaon process, both of which 
are more amenable to exact numerical evaluation. 

The (j)'s are computed recursively using the relation from stan- 
dard Markov chain theory: 


( 2 . 6 ) 


(})(0) = TT(0) 

(j)(n + 1) = (j)(n)P , n ^ 0 , 

(Note that Equation (2.6) involves only nonnegative numbers, a fact 
that contributes to numerical stability of the algorithm.) The matrix 
P is usually sparse and thus the above matrix multiplication should be 
performed by an appropriate algorithm. Such a multiplication algorithm 
is described by Gross and Miller [3]. The number of operations in this 
algorithm is proportional to the sum of the number of states and the 
number of transitions, e.g., 49 for the system of Figure 2. The programs 
in the appendices use this multiplication algorithm. 

In short, the standard randomization computational algorithm com- 
putes A and P from the generator Q using (2,1) and (2.2), respec- 
tively. It computes the truncation point T(0,t) from (2,4), then the 
cl)(n)'s using (2.6) recursively, accumulating in Equation (2.5) to give 
TT (t) . This algorithm was applied to ten versions of the model in Fig- 
ure 2. The results are summarized in Section 4. 


3, AN ACCELERATED ALGORITHM USING SELECTIVE RANDOMIZATION 

A close investigation of the standard randomization algorithm 
and the model of the three-component system in Figure 2 suggests a way 
to speed up the algorithm for this kind of model. In the three-component 
model states 0 , D , and DD have very long mean holding times because 
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the component Ejillure rate X is very small. All other nonabsorbing 
states have much shorter holding times. The absorbing states have infi- 
nite holding times. The process {X(t) , t ^ 0} spends most of its time 
in the states with long holding times. The Markov chain {Y^, n = 0,1,...} 
tends to sit in these states for many occurrences of the underlying Poisson 
process {N(t), t ^ 0} , making a null transition at each occurrence. 

By eliminating the computation involved in the’se null transitions for 
the states with the longest holding times, the speed of the algorithm 
can be Increased, 

Consider a modification of the model for the three-component system 
of Figure 2. The states AD , DD , ED , DD , FD , and DDA are each 
split into two states in order to distinguish whether or not the first 
fault is detected before the second fault occurs. The modified model is 
shown in Figure 4; it has 26 states and 42 transitions. This modification 
reveals (in Figure 5) a special block tree structure which can be ex- 
ploited in an accelerated randomization algorithm. The structure con- 
sists of the process alternating between states with long holding times 
(Si* S^, S^, and Sg in Figure 5) and short holding times (S 2 and S^) , 
not returning to any subset after leaving it, and finally being absorbed 
into a terminal set of states (Sg, S^, S^, . Larger, more realistic 

models of many fault-tolerant systems will tend to have this same struc- 
ture. Such a model is depicted in Figure 6. States with no undetected 
faults will have long holding times while those with undetected faults 
will be sliort. (Systems that contain processes with significantly dif- 
ferent time scales are called "stiff" in the literature on differential 
equations.) 
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Figure A Three~coinpoi\ent system reliability 
model: modified state space with 
tree structure. 


The accelerated randomisation algorithm is based on a semi-Markov 
process representation of the Markov process {x(t), t > 0} . (Ross 
[7] presents Markov processes as a special case of semi-Markov processes.) 
This representation is also an extension of an idea called "selective 
randomisation" by Mel.'unud and Yadin [6]. Soioctive randomisation is 
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DDA' 



S 


7 


S 


6 



Figure 5. — Block structure of modified state space of 
three-component system reliability model. 


similar to randomization except that the Markov chain is randomised by a 
Poisson process only wliilo it is in a subset of tlie state space 5 

The Markov chain is given arbitrary exponential holding times for tlic 
set S* = S - S|^ of exceptional states. (In the model of Figures 4 and 
5j Sj^ — S2 ^ Syj .) Let 
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^ System 



Failure 

—Lz 




All Faults 
Detected 



m 

^ System 

B 

Failure 


All Faults 
Detected 


Undetected 

Faults 


Figure 6. — Tret’~.iikc state space and transitions for a general model of 
fault-tolerant computers. (The initial state and states with 
all faults detected have very long mean holding times. States 
with undetected faults have very short mean holding times. 
System failure states have infinite holding times.) 


A - max q = max q , (3.1) 

sgSj^ ^ seS 

and define a substochastic transition matrix P* , 





A 


i e S„ 


>1i " ‘ 


" 'i.j 

q^ > 0, 1 E S'" 

(3.2) 

'.I 

'li 


0 

V 


q^ = 0. i £ S* . 
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Then {X(c), t > 0} can be characterized as a semi-Markov process with 
transition matrix P* and exponential holding times with rates {r^, 
i G S} , whore r^ « A for i ^ and r^ = for i e S* . (Thus 
the process is uniquely determined by tt( 0) , P* , A , and {q^, i G Ss'f} .) 
Let ^ “ 0,1,,.,} be the embedded Markov chain with substochastic 

transition matrix P* , and let N*(t) equal the number of transitions Of 
Z in [0,t] , noting that Z may make transitions i + i for i G Sj^ but 
not for i G S* . Denote the transient probability vectors for Z by 
i|{(n) , i.e., r|f(n) = (n) , , where (n) “ P(Z = s) . 

mL m o LI 

The processes Z and N* are dependent. The selective randomization 
formula is 


PCx(t) “ s3 = J pCx(t) = s, N*(t) » n} 
n=0 

CO 

“ I = s, N*(t) = n} 

n=0 

(3.3) 

CO 

= I P(Z = s) l’CN*(t) = n I Z = s3 
n=0 ” 


00 


= I K(n) pCn*(c) = 

n=0 ® 


n 


Z 

n 


s} 


The accelerated randomization algorithm is based on Equation (3.3). The 
]|i's can be comptifcd recursively, 


i|i(0) = TT(0) 

^(n + 1) = ^(n)P* 

In addition the quantltie.s PCN*(t) = n | Z^^ = 
the infinite series in Equation (3.3) imist be 
desired numerical accuracy. 


(3.4) 


.s3 must be computed, c* id 
truncated to achieve the 
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We note, for each subset of S in Figure 5, that 

PCN*(t) = n I = s3 , s e , takes a constant value; for 
1 = 1 , 2 , ... ,10 , 

^i “ PCN*(t) = n I Z^ = s} , s e , 

for n and s such that = s) 0 . Thus to use Equation (3,3) 

it suffices to compute A. (t) , i = 1,2,,.., 10 . As an example, con- 

JL I Ti 

sider A, (t) ; this also equals the probability of {N*(t) = n} given 

^ j 11 

{x(t) e s^} occurs. A typical sample path depicting this situation is 
sho™ in Figure 7 . In order to compute the probability of n occurrences 
in [0,t] we revert to the standard randomization construction; the 
holding time in state 0 has an exponential distribution with rate 2X . 

In the standard randomization, Pq 0 ~ ^ = (A - 3X)/A and 

Pq ^ “ 3X/A, and the transition to A will accur on th^ (i + l)st 
occurrence of the underlying Poisson process with probability 
(3X/A)d(A - 3X)/A3^ , i = 0,1,2,... . Similarly the holding time in 
state D has an exponential distribution with rate 2X and the process 
will leave D after being there for exactly j occurrences of the un- 
derlying Poisson process with probability (2X/A)(|(A - 2X)/A}'^ , 
j = 0,1,2,... . Consequently ut.are are many ways that {N*(t) = n) 
can occur, depending on the number of occurrences of {N(t), t > 0} 
that happen while {X(t), t > 0} is holding in 0 or D. Combining all 
these facts gives 
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Holding in Poisson Holding in Poisson 

State 0, Process State D, Process 

(rate A) T (rate A) 


Figure 7. — A typical sample path of process on modified 
state space for three-component system. 


■ Jo jo ^ ^ 

= r ? a lltnY M ( a-2X )J 

j~0 ^ ^ * A \ A / (n+i-Pj)! 


6A (f A r -2Xt -At [(A-2A)t]'^1 


- (^)"" [■ 


^ ^ -3At -At [(A-3.\)t]'^“ 
e - e A k! _ 


k=0 


T!icn analysis of A, (t) for other sets is based on the same principle 

1 1 n 

These quantities are: 
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“ T 4 '3'"> ’ " ^ ’ 

A 3 = 3a^“^ f^Cn) - 3a""^ f 3 (n) , n > 2 , 




- 2 X . 

" A ^3, 

n<'^) • 

n ^ 

2 

i 




„<o 

. ri 

■3 

= ^^1 

q(n) - 

6a;-2 

^2 

M 

+“ 333 "^ 

13 ( 0 ) 1 n > 3 . 


f “ 







(3.5) 

^6 

n<t) 

»n 

= 

- n -2 
- 3a^ 

« 2 <-) 

+ 

^4' 

£ 3 ( 0 ) 

> n > 3 > 



= fo(n) 

0 n-3 
- 3aj^ 

f j^(n) 

+ 

34 

f 2 (n) 

- 33 ^ f 3 (a) > n > - 

^8 

(t) 

,n' 

= fo<n) 

n -1 
- ^3 

£ 3 ( 0 ) 

> 

n 

> 2 , 



*S,n<« - f *l'"> - I ^ ^ ■ 

"lO.nC^) “ ^0<"> - f 'l<"> + i -r' ^ ^ ■ 

where, for 1=0, 1,2, 3, 

^ , -iXt -Xu r(A-iX)ul’" 

f (n) = e - z e 

^ k .=0 

A 

^1 “ A-iX ' 

These equacions provide Uhe required probabilities for Equation (3.3). 

To complete Uhe specification of the accelerated randomization 
algorithm it is necessary to give a truncation rule for the iniinite 
series in Equation (3.3). Note that 
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A 10 

pCn*(u) < Q ^ I I l>Cx(l:) G S., N*(U) = tO 

n«0 i«l ^ 


I 10 

I 

n«0 


X/ XU 

I ■’«„ V 

«0 i<=l ’ 


(3,6) 


a 10 

I I h I 

n=0 i«l seS^ 


and define 


T* (£,!:) = min 


A: PCN*(t) < si) > 1 - z 


(3.7) 


Thus, to summarize the accelerated randomization algorithm: the i/>'s 

are computed recursively using Equation (3.4) and the A's are computed 
using (3.5) with the products being accumulated in Equations (3.3) and 
(3.6) until the truncation point T* of Equation (3.7) is met, at which 
point the algorithm terminates, yielding probabilities which are accurate 
to within e of the exact values. 

The accelerated randomization algorithm was programmed (see the 
appendix for FORTRAN listing) for the model of Figure 4 and executed for 
the same ten versions of the system as the standard algorithm. Results 
are summarized in Section 4. 


4. COMPUTATIONAL RESULTS 

The standard randomization algorithm and the accelerated randomi- 
zation algorit)mi liave been programmed in FORTRAN for the system depicted 
in Figures 2 and 4, respectively. Listings of the programs appear in 
the appendix. Ten cases were run with different sets of parameter values. 
The input values for the different cases are given in Table 1. (Note 
that cases 3 and 4 and cases 5 and 6 arc identical except for the user 
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CASES COMPUTED 


Case 




Parameters 





ot 

P 

6 

Y 

p 

q 

Si 

E 

1 

lo"^ 

10 

1 

10 

100 

100 

.99 

1 

10"^ 

2 

lo"^ 

10 

1 

10 

100 

..100 

.99 

10 

lO"^ 


-3 








-9 

3 

10 

10 

10 

10 

100 

100 

.99 

1 

10 ^ 

4 

_3 

10 

10 

10 

10 

100 

100 

.99 

1 

10"^ 

5 

lO”^ 

10 

10 

10 

100 

100 

.99 

10 

10”® 

6 

10-^ 

10 

10 

10 

100 

100 

.99 

10 

10-’ 

7 

lo"^ 

100 

100 

100 

10* 

10^ 

.99 

1 

10”® 

8 

10"^ 

100 

100 

100 

10^ 

lo"^ 

.99 

10 

10”® 

9 

io“^ 

10^ 

lo'^ 

10^ 

lo'^' 

10^ 

.99 

1 

10"® 

10 

lO"^ 

100 

100 

100 

10^ 

10^* 

.99 

1 

io“® 

specified 

error 

bound, 

e .) 

The programs were run 

on The George Wash- 


ington University's IBM 370/4341. 

The execution times oO these randoinizafiioTi programs are propor- 
tional to the product of the truncation points (T or T*) and the sum of 
the number of states plus the number of transitions. The accelerated 
version requires mor‘e CPU time for each term in the randomizncion Eor- 
mula (3.3) because the weights A. (t) , i = 1,2 10 require more 

•L ) IT 

computation time. However, for systems \-?ith large state spaces this 
will be insignificant compared to the calculation in Equation (3.4). 
Tlius performance is more accurately predicted by the number of terms 
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multiplied by 18 + 31 = Cor the standard algorithm and 26 + 42 «= 68 
for the accelerated algorithm. The number of terms (truncation point) 
and CPU times in seconds are summarized in Table II. In most cases the 
accelerated algorithm appears to be far superior. 

The actual probabilities computed are presented in Table III, 

The probabilities listed are for the accelerated modification. The 
probabilities from the standard algorithm agree completely with these 
numbers and may be recovered from Table III by summing the probabilities 
for the split states, e.g,, P(AD) + P(AD') , 


TABLE II 

PERFORMANCE OF RANDOMIZATION ALGORITHMS: 
NUMBER OF TERMS REQUIRED AND CPU TIME 


Case 

Standard Algorithm 

Accelerated Algorithm 

No . Terms 

CPU Seconds 

No . Terms 

CPU Seconds 

1 

194 

3.20 

159 

6.50 

2 

1424 

19.54 

1233 

48.81 

3 

204 

3.12 

130 

5.40 

4 

193 

3.06 

95 

4.05 

5 

1522 

20.02 

206 

8.39 

6 

1492 

20.23 

143 

5.79 

7 

1522 

20.01 

174 

7.16 

8 

13656 

174.51 

208 

8.40 

9 

4383 

56.43 

95 

4.03 

10 

1522 

19.59 

143 

5.86 



STATE PROBABILITIES C0^^PUTED 
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5. CONCLUSIONS 


ORiGlNAL PAGEjS 


Randomization appears to be a good way to compute reliabilities 
for Markovian fault- tolerant computing systems with state spaces of mod- 
erate size. Gross and Miller [4] have solved Markov processes with 
20,000 states and 200,000 transitions using the standard randomization 
procedure. It is certainly feasible to use the approach on Markovian 
models of fault-tolerant systems of comparable* or even larger size. 

The accelerated randomization algorithm gives a significant sav- 
ings in CPU time for most examples. There should be an even greater 
savings for larger systems. Furthermore, this accelerated Implementation 
is applicable to any passage time problem, the exceptional set S* being 
the target states (with holding times set to infinity). This has appli- 
cation in computing fault-recovery-time distributions for fault-tolerant 
systems . 

The randomization algorithm is quite easy to implement. The main 
difficulty encountered in larger systems would be generation of the Q 
matrix. It is necessary to have an automatic way for the computer to 
genejrate Q or a sparse representation of it. Fortunately, the SERT 
methodology (see Gross and Miller [3]) can be applied to models of 
fault-tolerant systems to overcome this difficulty. 

The usual approach to computing transient probabilities for 
Markov processes is solution of the Kolmogorov forward equation 

tt' (t) = Tr(t)Q , t > 0 . 

This is an initial value system with Tr(0) given. There are two general 
approaches; (i) mimcricnl integration techniques such as Uunge-Kutta, 
predictor-corrector, etc., and (ii) exponentiation t'Jr(t) = 7r(0)e^*'] by 
computing the spectrum, computing the Taylor series, or other means. 
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The randomization technique has a distinct advantage over these ap- 
proaches in that a bound on the global error can be set by the user, 
and it is achieved with certainty. Furthermore, Grassmann [2] has shown 
randomization to be more efficient for some queuing systems. 

Another advantage of the randomization approach is that it is a 
"computational probability" technique. Computational probability is an 
emerging discipline concerned with numerical s'olution of applied prob- 
ability problems. The probabilistic structure of the model is exploited 
to obtain efficient numerical algorithms and to evaluate the performance 
of algorithms. In this particular application, probabilistic reasoning 
led to the accelerated algorithm. Another benefit of the probabilistic 
analysis is that Equation (2.3) just involves nonnegative numbers, 
creating numerical stability. Finally, the probabilistic point of 
view leads to efficient numerical algorithms for computing other quanti- 
ties of interest, for example, occupancy time distributions and expecta- 
tions which can be used in a performability analysis. 

Over all, it appears that tlie randomization technique is a very 
promising methodology for calculating reliabilities and related quan- 
tities for ^farkovian fault-tolerant computing systems. 
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THIS PROGRAM IS SPI-CIALIZED FDR A PARTICULAR 
HARKOV PROCESS UITH : 

18 STATES , 32 TRANSITIONS , 6 ABSORBING STATES , 

AND ONE INITIAL STATE ( E.D. STATE N0.1 ) 

THE ABSORBING STATES ARE : STATES NO- 6,7,8,13. 7,10 


INTEGER TSTAR<50) 

DOUBLE PRECISION RSTAR <20 ) , PSTAR ( SO ) 

DOUBLE PRECISION PH I OLD ( 2S ) , PH INI* U ( 2S ) ,PR0(2'5) ,RLOPRO<25) 
DOUBLE PRECISION RLAMDA, ALPHA , BETA , DELTA , RHO , OUE , EPS! LO 
DOUBLE PRECISION ERROR , TO , LT , TERH , TRSUH , PRLT , RLAHAX , RK , COUN 


READIN INPUT PARAMETERS 

READ, RLAHBA, ALPHA, BETA,BELTA,RHO,QUE,EPSILO 


PRINT, 


LAHDA 

•',RLAMBA 

PRINT, 

✓ 

ALPHA 

% ALPHA 

PRINT, 


BETA 

,beta 

PRINT, 


DELTA 

■',BEI_TA 

PRINT, 


RHO 

' , RHO 

PRINT, 

y 

OUE 

" ,oue 

PRINT, 

y 

EPSILON: ^ , EPSIL(D 

READ, 

TO 



READ, 

ERROR 


PRINT, 

y 

ERROR 

TERROR 


READ (5,7) (TSTAR( I) , 1=1 ,50) 
7 FORMAT ( 1 51^/1514/1514/514) 

NN = 1 8 


CONSTRUCT THE RSTAR VECTOR 


RSTAR (1) = 

“ 3 =♦: RLAMDA 


RSTAR (2) = 

3 ■* RLAMDA 


RSTAR ( 3 ) 

- ^ ALPHA- • RHO 

•+ -DELTA- +”2 

RSTAR<4) = 

ALPHA 


RSTAR(5) = 

RHO 


RSTAR(6) = 

DELTA 


RSTAR(7) = 

2 ^ RLAMDA 


RSTAR(S) = 

- ( DETA ••• 2 ii: 

RLAHBA ) 

RSTAR(9) - 

BETA 


RSTARM 0) 

= 2 =1= RLAMDA 


RSTAR ( 1 1 ) 

= - < EPSILO + 2 

■i: RLAHBA ) 

RSTAR ( 12) 

= OUE =t- EPSILO 


RSTAR ( 13) 

{ 1.0 - OUE ) 

=»= EPSILO 

RSTAR < 1 4) 

== 2 RLAMDA 


RSTARdS) 

-• - 2 =»• RLAMDA 


RSTAR ( 16) 

= 2 -1= RLAMDA 


RSTARM 7 ) 

-0.0 


RSTAR ( 1 8 ) 

= 0.0 


RSTAR ( 1 9) 

= 0.0 



=f:--RLAMDA-)^ 
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^ ORIGINAL PAGE KS RSTAR(20) = ~ ( Bl-Tfl ALPHA -h RHO + DELTA + RLAMRA ) 
OF POOR QUALITY RSTAR(21) = BETA 

RSTAR<22) = alpha 
RSTAR( 23) a RHO 
RSTAR<2^) = DELTA 
RSTAR(2‘o) « RLAMDA 

RSTAR(26) t= - ( 2 fJ; BETA + RI..AMBA ) 

RSTAR(27) == 2 * BETA 
RSTAR(28) = RLAHDA 

RSTAR(29) - - < BETA* EPSILO + RLANBA ) 

RSTAR(oO) = BETA 

RSTAROI) = OUE T EPSILO 

RSTARCJ2) « < 1.0 QUE ) >|i EPSILO 

RSTARC33) = RLAHDA 

RSTAR<3^) = - < BETA ••• RLAHBA ) 

RSTARC35) = BETA 
RSTAR(36) = RLAMBA 
RSTAR<37) -- 0.0 

RSTAR(38) = - ( ALPHA + RHO + DELTA + RLAMBA ) 

R5TAR(39) = ALPHA 

RSTAR('IO) « RHO 

RSTAR(AI) DELTA 

RSTAR<^2) = RLAMBA 

R3TAR(43) = - ( EPSILO + RLAMBA ) 

RSTAR<A‘I) = OUE T EPSILO 
RSTAR(45) - (1.0 - CIUE > EPSILO 
RSTAR<46) = RLAMDA 
RSTAR(47) = - RLAMBA 
RSTAR(48> s RLAMBA 
RSTAR(49) ••= 0.0 
RSTAR(SO) =0.0 
C 
0 

BO 18 1*.-1,NN 

— 18 PHIDLiiCn •= 0-0 

PHI0UK1) •" 1.0 

0 

C COMPUTE RLAMfVX 

C 

RLAMAX == RSTAR( 1 ) 

BO 19 1=2,50 

IF ( RSTARCn .LE. RLAMAX ) RLAMAX = RSTAR(I) 

19 CONTINUE 

RLAMAX = RLAMAK 
PRINT, •' RLAHAX : '■,RLANAX 
C 

COIJM = i;iLQC< RLAMAX n: TO ) 


23 


C 
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C CONSTRUCTION OF THE PSTAR UECTOR 

C 

no 20 1=1,50 

20 PSTAR(I) = RSTAR(I)/RLAHAX 
DO 21 1=1 ,50 

IF ( PSTAR<n .GT. 0.0 ) GO TO 21 
PSTAR(I) = 1.0 + PSTARC) 

21 CONTINUE 
C 

C 

c 

N = 0 

LT = - RI.AMAX -f- TO * s 

IF < LT .GT. -120 ) THEN DO 
TERM = nr£XP< LT ) 

ELSE DO 
TERM = 0.0 
END IF 

TRSIJM = TERM 
DO 27 1=1, NN 

27 PRO(I) = PHIOLIKI) * TERM 

IF ( TRSUM .GT. < 1,0 - ERROR ) ) GO TO 100 
29 CONTINUE 
PRLT = LT 
K = K + 1 
RK = K 

CALL TRpNC<COUN, RK, PRLT, LT, TERM) 

CALL EUAL(NH,PHIOLli,PGTAR,TSTAR,PHINEU) 

DO 33 I=1,NN 
PHIOLDTI) = PHINEU(I) 

33 PRO(I) = PRO a) + PHINEU(I) TERM 
TRSUM = TRSUM + TERM 

IF < TRSUM .GE. ( 1.0 - ERROR ) ) GO TO 100 
GO TO 29 
100 CONTINUE 

DO 110 1=1 ,NN 

RLOPRO(I) = DL0G1 0(F'ROn ) ) 

--tlO CONTINUE 

URXTE(6,130) 

130 FORMAT (••• 'rSX,'' TSTAR ■',20X,^ RSTAR ^SOX," PSTAR ',//) 

DO 135 1=1,50 

URITE(6,1 ^0) TSTARU ) ,RSTAR<I) ,PE-iTAR(I ) 

HO FORHATt-' , 7X , IN , 7X ,F28 . 23 , 6X , F20 . 23 ) 

135 CONTINUE 

UR1TE(A,150) K 

150 FORMATS' •'//■' H " NUMBER OF TERMS SUMMED : -’,1 1 7 ) 

WRITE (6, 160) TO 

160 FORMATC-' '(/' ' TIME OF INTEREST H,F10.6> 

URITE(6,162) 

162 FORMAT '////' ' THE LOGARITHM OF THE PROBABILITIES IS : 

DO 165 1=1, MN 

URITE(6,170) I,RLUPRO<I) ,PRO(I) 
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170 FORMAK'' '■//'' '■jSX,’' LOGPRQ ( " , 13 , 
165 CONTINUE 
STOP 
END 
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^,F30.26,5X,F30.26) 


SUBROUTINE EV AL < NN , PHIOLD , PSTAR , TST AR, PHI MEU ) 
DOUBLE PRECISION PHIOLD ( 25 ) , PSTAR < 50 ) , PHI NEU ( 25 ) 
INTEGER TBTAR(50) 

DOUBLE PRECISION PHIJ 
BO 33 J=1yNN 
PHINEU(J) =0.0 
33 CONTINUE 
I = 0 

DO 2 J'-l ,NN 
I = I 1 

PHIJ = phiol:o(J) 

PHINEU(J) = PHIJ f».' PSTAR(I) + PHINEU(J) 

MJ = TSTARd) 

IF ( NJ .EQ. 0 ) GO TO 2 
no 1 K=1,MJ 
I = I 1 
LJK « TSTAR<I) 

1 PHINEVMLJK) = PHINEU(LJK) + PHIJ PSTARd) 

2 CONTINUE 
RETURN 
END 


SUBROUTINE TRUNC ( COUN , RK , PRLT , LT . TERN ) 
DOUBLE PRECISION CDUN , RK , PRLT , LT , TERM 
LT = PRLT + CDUN - DLOG( RK ) 

IF < LT .GT. “120 ) THEN DO 
TERM = DEXP( LT ) 

Etr&E' Dii 

TERM =0.0 
END IF 
RETURN 
END 
C 

C END OF PROGRAM 1 LK 
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C THIS PROGRAM IB SPECIALIZEP FOR A PARTICULAR 

C MARKOV PROCESS UITH : 

C 26 STATES , A2 TRANSITIONS , 10 ASSORBING STAT) 

C 

IMPLICIT REALtS(A-H, O-Z), INTEGER<I-N) 

INTEGER TSTAR<68) 
niMENSION RSTAR(68) ,PSTAR(A£)) 

DIMENSION PHIOLIKIiA) ,PH;i:NI-IJ< 24) ,PR0(2A) ,RLOPRO(2<i) 
DIMENSION PIK(1 0) ,PGI< 1 0) ,DLA(4) 

DIMENSION COUN( 4) ,ST(4) ,PRST( 4) ^PRLT( 4) 

DOUBLE PRECISION LT(4) 

DOUBLE PRECISION ALG ( 3 ) , FI ( 4 ) 

DATA XL/1.00D-A2/ 


C 

C 

C 


C 

C 

C 


C 


Q 

C 

C 


READIN INPUT PARAMETERS 

READ, RLAMDA, ALPHA, BETA, DELTA, RHO,QUE, EPSXLO 

PRINT, ^ LAMDA :'',RLAMDA 

PRINT, ALPHA i'*, ALPHA 

PRINT, •' BETA :-',BETA 

PRINT, •' DELTA s^J^ELTA 

PRINT, RHO !’'.RHQ 

PRINT, ■' RUE sSQlJE 

PRINT, ■' EPSILONs , EPSILO 

READ, TO 

READ, ERROR 

PRINT, ERROR ERROR 
READ, NN 

READ (5,1 ) (TSTARCn ,I“-1 ,68) 

FORMAT (50 1 5/50 i5/wI5) 

THE riLA PARAMETER VECTOR 


DLA( 1 > 

= 3 RLAMDA 


DLA(2) 

« 2 RLANBA 


DLA(3) 

= RLAMDA 


DLA(4) 

= 0.0 


DXP1 

»EXP( ( riLAd) 

>1= TO) ) 

PRINT, 

■” DXfM :'',DXPI 


DXP2 - 

DEXP( - ( DLA(2) 

!l= TO)) 

PRINT, 

■' DXP2 :',riXP2 


DXP3 = 

DEXP( - ( DLACJ) 

TO ) ) 

PRINT, 

•' DXP3 :‘',DXP3 



CONSTRUCT THE RSTAR OECTOR 


RSTAR< 1 ) 
RSTARO) 
RSFARCn 
RSI ART 4 ) 


0.0 

0.0 

- ( ALPHA RHU tiELTA + 3 RLAMDA ) 
ALPHA 
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RSTAR(5) = RHO 

RSTARC6) « DELTA 

RSTAR(F) = 2 RLAMDA 

RSTAR(8) « “ < BETA 2 '<= RLAMDA > 

RSTAR(9) = beta 


RBTAR( 1 0) 


2 •)' RLAHBA 

RSTARd 1 ) 

ts 

- ( EPBILO + 2 RLAHDA ) 

RSTAR( 12) 

s:; 

aUE EPSILO 

RSTARd 3) 

tz 

( 1.0 - QUIi ) .1; EPSILD 

RSTARd 4) 

= 

2 =»! RLAfiDA 

RSTARd 5) 


-> ( ALPHA + BETA + RHO + DELTA + 

RSTARd 6) 

is 

ALPHA 

RSTARd 7) 


RHO 

RSTARd 8) 

- 

DELTA 

RSTARd 9) 

=5 

BETA 

RBTAR<20) 

rs 

RLAMBA 

RSTAR(21 ) 

ss 

- ( 2 t. BETA RLAMBA ) 

RSTAR(22) 


2 -I! beta 

RSTAR(23) 

!= 

RLAHDA 

RSTARC24) 

:= 

- ( BETA + ERSTLO + RLAMBA ) 

RSTAR( 25) 

:= 

UUE =<•• EPSILO 

RSTAR ( 26 ) 


BETA 

RSTAR(27) 

=r 

( 1.0 - QUE ) EPSILO 

RSTAR(2S) 


RLAHDA 

RSTAR(29) 

3 

- ( BETA + RLAMBA ) 

RSTAR<30) 


BETA 

RSTAR<31 ) 

s= 

RLAMBA 

RSTAR(32) 

3 

- ( ALPHA RHO BELTA + RLAMDA 

RSTAR(33> 


ALPHA 

"RSTAR<34) 

■*23. ' 

RHO 

RSTAR(35) 

u: 

RLAMDA 

RSTAR<36) 

=3 

DELTA 

RBTAR(37) 

;s 

- ( EPSILO + RLAMDA ) 

RBTAR(38) 

= 

( 1.0 “ OUE ) •(= EPSILO 

RSTAR v39) 

3 

RLAHDA 

RSTAR('IO) 

3 

CHJE •■(i EPSILO 

RSTAR< 41 ) 

i: 

0.0 

RSTAR < 42) 

::: 

0.0 

RSTAR < 43 ) 

is: 

- ( ALPHA -i- RHO + BELTA + RLAMDA 

RSTAR (44) 

3 

ALPHA 

RSTAR (45) 

3 

RHO 

RSTAR (46) 

;s 

DELTA 

RSTAR ( 47 ) 

£2 

RLAMBA 

RSTAR ( 48) 

'.n 

■" ( BETA + RLAMDA > 

RSTAR (49) 

S3 

BETA 

RSTAR (50) 

3 

RLAMBA 

RSTAR (51 ) 

:= 

- ( EPSILO + RLAMDA ) 

RSTAR (52) 

33 

CUIE Hi EPSILO 

RSTAR (53) 

S3 

( 1.0 - nUE ) ip EPSILO 

RSTAR (54 ) 

S3 

RLAH BA 

RSTAR (55) 

U 

0.0 

RSTAR (5ft) 

53 

0.0 


tLAMDA ) 
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RSTAR(57) 

tz. 

0.0 

RSTAR?'58) 

til 

0.0 

RSTAR(59) 

s 

0.0 

RSTAR( 60) 

til 

0.0 

RSTAR(61 ) 

:a 

0.0 

RSTAR(62) 

:= 

0.0 

RSTARC63) 


0.0 

RSTAR< 6^) 

;r. 

0.0 

RSTAR(65) 

IS 

0.0 

RSTAR<66) 

=i 

0.0 

RSTAR(67) 


0.0 

RSTAR(68) 

tz 

0.0 
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c 

c 

c 


1 1 


c 

C-. 

c 


1 y 


21 


c 


DO S 1=1,26 
PHIDLXra) =: 0.0 
CONTINUE 
PHIQI.IK 1 ) = 1 .0 

COMPUTE RLAHAX 

RLAMAX = RSTAR<1) 

DO 11 1=2,68 

IF ( RSTARd) .LE. RLAMAX ) RLAMAX « RSTARU’.) 
CONTINUE 

RLAHAX = ~ RLAHAX 
PRINT, •' RLAHAX :’',RLAMAX 


C 0 H S T K ij C T i ( j H 0 i" ‘ T H i-i P 8 T A k C T u K 


DO 17 1=1 ,68 

PSTAR(I) = RSTARd )/RLAHAX 
DO 21 1=1,68 

IF ( PSTARd) -GT. 0.0 ) GO TO 21 
PSTARd) = 1.0 + PSTARd) 

CONTINUE 
PSTARd) = 0.0 
PSTAR(41 ) =0,0 
PSTAR(55) = 0.0 
PSTART57) =0.0 
PSTAR(58) = 0.0 
PSTARd'?) =0.0 
PSTAR(60) =0.0 
PSTARd>1) = 0.0 
PSTAR(62) =0.0 
PBTAR(63) •= 0.0 
PSTAR(6'*) = 0.0 
PSI'AR<63) = 0.0 
PSTAR(66) 0.0 
PSTAR(68) = 0.0 
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ALG(I) = RLAMAX / ( RLAMAX - BLA(1) ) 

ALG(2) « RLAMAX / < RLAMAX - HLA(2) ) 

ALG(3? =* RLAMAX / ( RLAMAX - »LA<3) ) 

C 

C START CALCULATIOMS X = 0 .... 

C 

K = 0 

TRSUM = 0.0 

DG 24 I-l ,26 s 

PRG<I> = 0.0 

24 CONTINUE 

PROM ) = riEXP( - < ULAM ) =i! TO ) ) 

TRSUM = PROM) 

IF ( TRSUM .G£. ( 1.0 •“ ERROR )) GO TO 333 

C 

DO 26 1=1 ,4 
ST<I) = 0.0 
PRSTM) = 0.0 
LT<I> = 0.0 
PRLTM) =0.0 
FKI) = 0.0 

26 COUN<I) = HLOG< (RLAMAX - DLA(I)) •+ TO) 

C 

K = 1 

25 CONTINUE 

DO 27 1=1,10 

27 PIK(I) =0.0 

CALL EUAL(NN,PHIOLD,PSTAR,TSTAR.PI-!INEU) 

-e 

C CALCULATE THE STM!, TO, K) QUANTITIES 

C 

IF ( K .GT. 1 ) GO TO 31 

DO 29 1=1 ,4 

LT(I) = " RLAMAX * TO 

IF ( Lf(I> „(3T. -120 ) THEN DO 

ST(I) = DEXP{ LT(I) ) 

ELSE DO 
STM) = 0.0 
END IF 

29 CO NT IN LIE- 

PRINT, ■' STM,TO,t) :",STM) 

PRINT, ■' ST(2, TO, 1) :',ST(2) 

PRINT, ST(3,T0,1) !-',ST<3) 

PRINT, ■' ST(4,T0,n :-',ST(4) 

GO TO 33 
31 CONTINUE 
RK = K 
1 = 1 

C A 1. 1„ 3 LI M ( I , 1^ K , P S -1- , R L T , C 0 U N , T E R N , L. T , S 'I' ) 

1=2 

RK = RK +1.0 

C A L I. S tJ M < I , R K , P R S T , PR 1.T , C OU N , TE RM , 1. T , ST ) 
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I => 3 

RK = RK + 1 .0 

CALL SUM<I,RK,PRSTrPI\’LT,COU»,TERMrLT,ST) 
I = 

RK 5 RK + 1.0 

CALL SUMa,Rf{,PRST,PRLT,COUN,T£RH,LT,ST) 
33 CONTINUE 
C 

Fid) = nxpl - STd ) 

FI (2) a DXP2 ^ STC2) 

FI (3) = DXP3 ~ ST(3) 

FI (*U a 1.0 - ST(V|) 

(3 CALCULATE THE PIK ^ S ‘1 

PIK<1) a 0.0 

CALL TlJOaiLA,RLAMAX,ALG,K,FI,XXX) 

PIK<2) a XXX 

CALL THREE <ALG,I(, FI, XXX) 

PIK(3) a XXX 

CALL FQIIR(DLA,RLAMAX,PIK,XXX) 

PIK(‘1) = XXX 

CALL FIVE(ALG,K,FI,XXX) 

PIk(5) a XXX 

CALL SIX(FI ,ALG,K,XXX) 

PIt((6) = XXX 

CALL S£OEN<FI ,ALG,K,XXX) 

PIK(?) a XXX 

CALL E I GH T ( F T , A LG , K , XXX ) 

PiK(B) a XXX 

CALL NINE(ALG,K,PI,XXX) 

PI Id 9) a XXX 

CALL TEN(I"I , ALG, K,XXXJ 

PIHdO) = XXX 


C EVALUATE THE CURRENT PRORAJULITIES 


C 

DO 3-1 1=1,10 

IF ( Plldl) .LE. XL ) THEN DO 
Plldl) a 0.0 
ELSE DO 

PIK(t) a PIKM) 

END IF 

3-1 CONTINUE 

DO 2U :d1,i'{N 

IF ( PIIINEUd) -LE. XL ) THEN na 
PHINEUd) = 0.0 
ELSE IH) 

PHINEU(I) = PH 1NEU< I) 

END IF 

33 COimmjE 
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PROd) = PROM) ••• 
DO 30 

PROd) = PRO(I) + 
38 CONTINUE 

PROM 1 ) = PROM 1 ) 
DO 37 I = 

PROd) = PROd) + 
37 CONTINUE 

PR0<15) = PROMS) 
DO ‘10 1 = 16,1? 
PRO<I) ••= PROd) + 

40 CONTINUE 

PROM 8) = PR0<1B) 
DO ‘11 1 = 17,2-1 
PRO<I) = FROM) + 

41 CONTINUE 

PRO (23) = PRO (213) 
PRO (26) = PRO (26) 
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RHXNEUM) PIKM) 

PHINliUd) 'I' PI|{(2) 

+ PMINEWM 1 ) RIK(3) 

PHINEUd) 'i! P1K(4) 

+ PHINEUM5) •«* PIK(S) 

PHINEUd) -v- PII(<6) 

+ PHINEUM8) PIK(7) 

PHINEUd) >f« FIR (8) 

+ PHINEU(25) PII((?) 

+ PHrNEU(26) Hi PIKMO) 


CHECK IF YOU ACIEVED THE TRUNCATION POINT 


DO -1‘1 1 = 1,10 

44 PSId) = 0.0 
PSIM ) = PHINEUd ) 

DO 4S 1=2,10 

45 PSK2) = PSK2) + PHINEUd) 

PBI(3) = PHINEUMT) 

DO 46 1=12,14 

46 PSK4) = PSK4) + PHINEUd) 

PSK5) = PHINEUdS) 

DO 47 1=16,17 

47 PSK6) « P8K6) + PHINEUd) 

PSK7) = PHINEUM8) 

DO 48 1=19,24 

48 PSI(Q) = P8K8) 4- PHINEUd) 

P5K7) = PHINEU(2S) 

PSIMO) = PH1'NEU(26) 

C 

C 

DO 50 1 = 1,10 

50 TRSUH = TRSUH + PSI ( I ) -i* PIKd) 

IF ( TRyUft .GE. ( 1.0 - ERROR )) 00 TO 333 
K = |( + 1 
DO 51 1=1,4 
PRST(I) = SKI) 

PRLT(I) = LTd) 

51 CONTINUE 

DO 55 1 = 1 ,NN 

55 PHIOLD(I) = PHINEUd) 

GO TO 25 
333 CONTINUE 
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SUMHA ■= 0.0 

DO 73 1=1 ,NN 

SUMMA = SUMMA PRO (I) 

73 RLOPROd) = DL0G10(PR0a ) ) 

PRINT, ^ SUM OF PRODABILITII'S % SUMMA 
URJ.TE(6,130) 

130 FORMAT('' '//' TSTAR ’',20X,'' RSTAR SaOX," PSTAR ",//) 

DO 135 1=1,68 

URITEU,MO) TSTAR < I) ,R5TAR(I) ,PSTAR(I) 

HO FORHAK *' , 7X , H , 7X , F2B . 20 , 6X , F28 . 20 ) 

135 CONTINUE 

URITE(6,150) K 

150 FORMAT V/-' ' NUMBER OF TERMS SUMMED ;-',I7) 

URITE(6,160) TO 

160 F0RMAT1-' '//' ' TIME OF INTEREST :",F10.6) 

WRITE (6, 162) 

162 FORMAT(^ '■////'“ THE LOGARITHM OF THE PROBABU-ITIES IS s O 

DO 165 1=1, NN 

WRITE (6, 170) I,RLOPRO(I ) ,PRO(I) 

170 FORHAT<-' '//' ',5X,'‘ LOGPRO < , 13 H ) = •' ,1-30 .26 , 5X ,F30 .26 ) 

165 CONTINUE 
STOP 
END 
C 

C 

c 

-SUBROUTINE EVAL<MN,PHI0L,D,PSTaR;T5Tar,THINEW) 

DOUBLE PRECISION PHI0LD(26) ,PSTAR(68) ,PHINEW(26) 

INTEGER TSTAR(68) 

DOUBLE PRECISION PHIJ 
DATA XK/1 -OOE-62/ 

DO 77 J=i,NN 
PHINEW(J) = 0.0 
77 CONTINUE 
I =• 0 

DC) SI J = 1,NN 

I = I 1 

PHIJ = PHUILB(J) 

PHINEIKJ) = PHIJ ‘1= PSTAR(X) + PHINEIHJ) 

MJ = TSTAR(I) 

IF ( HJ .E«. 0) GO TO 81 
DO 79 KM=1 ,HJ 
1 = 1 + 1 
LJK = TSTAR (I) 

IF ( PHIJ .LE. XK .OR. PSTAR(I) .LE. XK ) THUN D(J 
PHINEUU-JK) = PHINEUaJK) 

ELSE DO 

PHINEU(LJK) = PIIINEWO.JK) + PHIJ + PSTAR(I) 

END IF 

79 CONTINUE 
G1 CONTINUE 
RETURN 
END 
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c 

SUBROUTINE TWO ( DLA , RLAHAX , AI.G , N , FI , XXX ) 

DOUBLE PRECISION DLA ( A ) , ALG ( 3 ) , F I ( 4 ) 

DOUBLE PRECISION RLANAX,XXX 
XXX = ALG(1)m>!|!K 

XXX = XXX DLA(I) Pin) / RLAHAX 
RETURN 
END 
C 

c 

c 

SUBROUTINE THREE ( ALG , K , FI , XXX ) 

DOUBLE PRECISION ALG ( 3 ) , F I ( -) ) 

DOUBLE PRECISION XXX 

XXX = < ALG (2 )=<=!(:< K ~ 1)) 3 Ht Fr<2) 

XXX = XXX - < ( ALG(1 )=!=■<=( K - 1 )) 'i- 3 FKl) > 

RETURN 

END 

C 

C 

c 

SUl^iROUTlNE FOUR(DLA,RLAMAX,PIK,XXX) 

DOUBLE PRECISION DLA ( *1) ,PIK ( 1 0 ) 

DOUBLE PRECISION RLAMAX,XXX 
XXX = DLA<2) H: PIN (3) / RLANAX 
RETURN 
END 


C 

SUBROUTINE FI VE < ALG , K , F I , XXX ) 

DOUBLE PRECISION ALG ( 3 ) , F I( 4 ) . 

DOUBLE PRECISION XXX 
IF ( K .GT. 1 ) (30 TO 1005 
XXX = 0.0 
GO TO 1006 

1005 XXX = ( AL0(3)i«-MK-2) ) -■i- 3 FI (3) 

XXX = XXX - ( ( ALG<2) -(-^(K-2) ) 6 FI(2) ) 

XXX = XXX + ( < ALG< 1 )*=MK-2) ) * 3 FK|) ) 

1006 CONTINUE 
RETURN 
END 

C 



C 

SUBROUTINE SI X (F T , ALD , K , XX X ) 

DOUBLE PrfECISlON FI ( 4 ) , ALG ( 3 ) 

DOUBLE PRECISION XXX 

IF ( K .GT. 1 ) GO TO 1011 
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XXX = 0.0 
RETURN 

toil XXX « [‘.T(4) “ ( { ALG<;>).tiMK-2) ) 3 |*I(2) ) 

XXX “ XXX + < ( Al.G( t )i*==MK-2) ) <»! 2 1-I(U ) 

RETURN 

ENti 

i: 

c 

c 

subroutine BEUEN(FI,ALG,K,XXX) 

DOUBLE PRECISION FI(n),ALG(3) 


DCIUBLE PRECl 

s 

ION 

XXX 



IF < K .GT. 

2 

) 

GO TO 1021 



XXX » 0.0 






RETURN 






1021 XXX = i-n'i) 


( 

( ALG(3) li'MK-3) 

) 

3 FI (3) ) 

XXX = XXX + 

< 

< 

AL(3(2)'i-‘MK-3) ) 

i|; 

3 F;t(2) ) 

XXX « XXX - 

( 

< 

ALG( n=*‘=M1(-3) ) 

!)! 

Fid) ) 


RETURN 

END 

C 

c - 

0 

SUBROUTINE El GUT ( F I , AUG , K . XXX > 
tiQUBLE •PkECISION FXi'i ) , ALiKS)* 

DOUBLE PRECISION XXX 

XXX - FI(A) - ( ( ALG< 1 )=inMK“1) ) FI(l) ) 

RETURN 

END 

C 

C 

SUBROUTINE NINE ( ALG , K , F I , XXX ) 

DOUBLE PRECISION ALG(3),PI<^) 

DOUBLE PRECISION XXX 

XXX " ( ( ALG(3)h.!MIC*-U ) +• 3 ms FI (3) ) / 2 
XXX = XXX - ( ( ALG( 3 FKU ) / 2 

RETURN 
END 
C 

c 

c 

SUBROUTINE TEH ( F I , ALG , K , XXX ) 

DOUBLE PRECISION F I ( <n , ALG ( 3 ) 

nolJDLF. PRECISION XXX 

IF C K .GT. I ) GO TO 1031 

XXX = 0.0 

RETURN 

1031 XXX - FI (‘1) - ( ( ALG< 3) ♦-MK-2) ) ts 3 m3) ) / 2 

XXX = XXX ( ( ALG( 1 FlU) ) / 2 

RETURN 
END 
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SUBROUTINE SUH(I,RK,PRST,PRLT,!:OUN,TERM,LT,ST) 

DOUBLE PRECISION PRST < 4 ) , PRLT ( 4 ) , COUN ( 4 ) , LT ( ^ , ST ( ^ ) 
DOUBLE PRECISION RK^TERM 
RK = RK - 1.0 

CALL TRMa , RK,PRLT,C0UN,LT,TERH) 
sra) « PRST(I) + TERN 
RETURN 
END 
C 

C 

C 

SUBROUTINE TRM (I , RK , PRLT , COUN , LT , TERM ) 

DOUBLE PRECISION PRLT < 4 ) , COUN < 4 ) , LT ( 4 ) 

DOUBLE PRECISION TERM,RK 

LTU) « PRLTCI) + COUN(I> - I.tL06(RK) 

IF ( LT(I) ,6T. - 50.0 > GO TD SOO 

TERN =0.0 

RETURN 

500 TERM = DEXP( LT(I> ) 

RETURN 

END 



